clear all;
clc;
tstart=1; %start time
tmax=200; %max time of simulation
time=tstart:tmax; %time vector

n_iterations=5;
%nparticles = 1000;
nparticles = 10;
npos=2*length(time)+3; %size of pos arrays = 2t+3
pos=1:npos;
nmid_rel=tmax+2; %relative mid position in the below a arrays which corresponds to zero position in the absolute sense
nstart_rel=0;%relative start position in the below a arrays which corresponds to least -ve position in the absolute sense
nend_rel=0;%relative end position in the below a arrays which corresponds to max +ve position in the absolute sense

a_n_u_prev=zeros(1,npos); %a's for spin up at each N at previous time
a_n_d_prev=zeros(1,npos);%a's for spin down at each N at previous time
a_n_u_cur=zeros(1,npos);%a's for spin up at each N at current time
a_n_d_cur=zeros(1,npos);%a's for spin down at each N at current time

%sigma=zeros(1,npos); %the standard deviation of the spatial distribution
sigma=zeros(1,length(time));
sigma_up=zeros(1,length(time));
sigma_down=zeros(1,length(time));

avg_pos=zeros(1,length(time));
avg_pos_up=zeros(1,length(time));
avg_pos_down=zeros(1,length(time));

sigma_count = 1;
% T=tmax/2:100:tmax;
T=time;
pos_counter=0;
a=0;
epsilon = 0.5;
p_n_tmax = zeros(1,npos);

p_n_u_tmax = zeros(1,npos);
p_n_d_tmax = zeros(1,npos);

p_n = zeros(1,npos);
p_n_u= zeros(1,npos);
p_n_d= zeros(1,npos);

sig=zeros(1,length(time));
sig_up=zeros(1,length(time));
sig_down=zeros(1,length(time));

%for all particles
% for p=1:nparticles
for iter=1:n_iterations
    display(iter);
    %form the initial a vectors for time = 0. consider all the particles to be
    %with spin up
%     a = 1/2;
    a = 1/(sqrt(2));
%     a = rand;
    a_n_u_prev(nmid_rel)=a;
    a_n_d_prev(nmid_rel)=(sqrt(1-a^2))*i;
%     a_n_d_prev(nmid_rel)=a*i;
    
    
    %for all times
    for t=time
%         display(t);
        pos_counter = 0;
        nstart_rel = nmid_rel - t;
        nend_rel = nmid_rel + t;
        if(t ~= 1)
%              a=1/sqrt(2);
             a=(1/(sqrt(2))) + epsilon*rand;
%             a = epsilon*rand;
        end;

        %for all possible positions for this particular time, fill up the a
        %arrays for the current time
        
        for pos_counter=nstart_rel:nend_rel

            a_n_u_cur(pos_counter)=  ( sqrt(1-a*a) * a_n_d_prev(pos_counter-1) )+ ( a * a_n_u_prev(pos_counter-1) );
            a_n_d_cur(pos_counter)=( sqrt(1-a*a ) * a_n_u_prev(pos_counter+1)) - ( a * a_n_d_prev(pos_counter+1) );

%             a_n_u_cur(pos_counter)=  (1/sqrt(2)) * (a_n_d_prev(pos_counter-1) + a_n_u_prev(pos_counter-1) );
%             a_n_d_cur(pos_counter)= (1/sqrt(2)) * (a_n_u_prev(pos_counter+1) - a_n_d_prev(pos_counter+1) );
            
            
            %record the probability densities of up,down and all particles
            p_n_u_tmax(pos_counter) = (abs(a_n_u_cur(pos_counter)))^2;
            p_n_d_tmax(pos_counter) = (abs(a_n_d_cur(pos_counter)))^2; 
            p_n_tmax(pos_counter) = (abs(a_n_u_cur(pos_counter)))^2 + (abs(a_n_d_cur(pos_counter)))^2; 
        end;
%         if(iter == n_iterations)

%             sigma(t) = sqrt( sum((pos-nmid_rel).^2 .* p_n_tmax) - (sum((pos-nmid_rel).*p_n_tmax)).^2 );
%             sigma_up(t) = sqrt( sum((pos-nmid_rel).^2 .* p_n_u_tmax) - (sum((pos-nmid_rel).*p_n_u_tmax)).^2 );
%             sigma_down(t) = sqrt( sum((pos-nmid_rel).^2 .* p_n_d_tmax) - (sum((pos-nmid_rel).*p_n_d_tmax)).^2 );
 
            sigma(t) = ( sum((pos-nmid_rel).^2 .* p_n_tmax) - (sum((pos-nmid_rel).*p_n_tmax)).^2 );
            sigma_up(t) = ( sum((pos-nmid_rel).^2 .* p_n_u_tmax) - (sum((pos-nmid_rel).*p_n_u_tmax)).^2 );
            sigma_down(t) = ( sum((pos-nmid_rel).^2 .* p_n_d_tmax) - (sum((pos-nmid_rel).*p_n_d_tmax)).^2 );
            
            
            avg_pos(t) = sum((pos-nmid_rel).* p_n_tmax);
            avg_pos_up(t) = sum((pos-nmid_rel).* p_n_u_tmax);
            avg_pos_down(t) = sum((pos-nmid_rel).* p_n_d_tmax);
%             sigma_count = sigma_count + 1;
%         end;
        
        %calculate the sigma 
        p_n = p_n + p_n_tmax;
        p_n_u = p_n_u + p_n_u_tmax;
        p_n_d = p_n_d + p_n_d_tmax;
        
        sig = sig + sigma;
        sig_up = sig_up + sigma_up;
        sig_down = sig_down + sigma_down;
        
        %set the a arrays of prev time with values of that of the current
        %time and a arrays of current time with zero values
%         if(t <tmax)
            a_n_u_prev = a_n_u_cur;
            a_n_d_prev = a_n_d_cur;
            a_n_u_cur = zeros(1,npos);
            a_n_d_cur = zeros(1,npos);
%         end;
        if(iter == n_iterations-1) 
            p_n_tmax = zeros(1,npos);
            p_n_u_tmax = zeros(1,npos);
            p_n_d_tmax = zeros(1,npos);
        end;
        
    end;
   
end;

p_n = p_n/(n_iterations*tmax);
p_n_u = p_n_u/(n_iterations*tmax);
p_n_d = p_n_d/(n_iterations*tmax);

sig = sig/(n_iterations*tmax);
sig_up = sig_up/(n_iterations*tmax);
sig_down = sig_down/(n_iterations*tmax);

% figure(1)
% plot(pos-nmid_rel,p_n);
% grid on;

% figure(2)
% plot(T,sigma,T,sigma_up,T,sigma_down);
% grid on;

figure(3)
plot(pos-nmid_rel,p_n,pos-nmid_rel,p_n_u);
xlabel('T:1000');
ylabel('P_t(n) for charge and up-spin');
title('P_t(n) for charge and up-spin, for \epsilon = 0.5');
grid on;

% figure(4)
% plot(pos-nmid_rel,p_n_tmax);
% grid on;

% figure(5)
% plot(pos-nmid_rel,p_n_u);
% grid on;


% figure(8)
% plot(T,avg_pos);
% grid on;
% 
% figure(9)
% plot(T,sig,T,sig_up,T,sig_down);
% xlabel('T:1000');
% ylabel('\sigma^2(t) for charge and up-spin,down-spin');
% title('\sigma^2(t) for charge and up-spin,down-spin, for \epsilon = 0.3 and a=1/sqrt(5)');
% grid on;
% 
% dlmwrite('sig2_0_8.dat',sig);
% dlmwrite('sig2_up_0_8.dat',sig_up);
% dlmwrite('sig2_down_0_8.dat',sig_down);